%----------------------------------------------
% Formatting and file path
%----------------------------------------------

clear;
clc;
format bank; % Set formatting to display numbers without e format

%----------------------------------------------
% Read in data
%----------------------------------------------
Folder = cd; 
Folder = fullfile(Folder, '../../Data/structural/');
cd(Folder); 
data = csvread('intermediate/data_for_matlab_stups.csv');
hhid=data(:,1);
l_bl=data(:,2);
h_bl=data(:,3);
tot_hrs_bl=data(:,4);
k_bl_post_trans=data(:,5);
hired_l_bl=data(:,6);
w_bl=data(:,7);
case_bl=data(:,8);
l_w2=data(:,9);
h_w2=data(:,10);
tot_hrs_w2=data(:,11);
k_w2=data(:,12);
hired_l_w2=data(:,13);
w_w2=data(:,14);
case_w2=data(:,15);
k_bl_pre_trans=data(:,16);
H=data(:,17);

N=length(hhid);    % Number of HHs
R=3650;            % Time endowment constraint
Hire_in_max=1400;  % Maximum number of hours that can be hired in

a=-1.88e-10;       % Estimated parameters
b=0.0009368;
beta=0.3499541;

uns_threshold=10380;

phi = 0.7;          % Ratio of hired in labour to market wage
rho = 0;            % Share of capital profit from liquidating


%----------------------------------------------
% Simplify terms
%----------------------------------------------

fk_w2 = a.*k_w2.^2 + b.*k_w2;
fk_bl = a.*k_bl_pre_trans.^2 + b.*k_bl_pre_trans;

%----------------------------------------------
% Calibrate psi and A parameters
%----------------------------------------------

psi_h = zeros(N,1);
psi_l = zeros(N,1);
A = zeros(N,1);

for i = 1:N
    % Mixing with hired-in labour in Year 2
    if case_w2(i) == 1
        A(i) = phi * w_w2(i) / (fk_w2(i)*beta*(l_w2(i) + hired_l_w2(i))^(beta-1));
        psi_h(i) = w_w2(i) / (phi*l_w2(i) + h_w2(i));
        psi_l(i) = phi^2 * psi_h(i);
    end

    % Livestock rearing only with hired-in labour in Year 2
    if case_w2(i) == 3
        if case_bl(i) == 6
            psi_h(i) = w_bl(i)/h_bl(i);
            A(i) = phi * w_w2(i) / (fk_w2(i)*beta*(l_w2(i) + hired_l_w2(i))^(beta-1));
            psi_l(i) = phi * w_w2(i) / l_w2(i);
        elseif case_bl(i) == 7
            A(i) = phi * w_w2(i) / (fk_w2(i)*beta*(l_w2(i) + hired_l_w2(i))^(beta-1));
            psi_l(i) = phi * w_w2(i) / l_w2(i);
        end
    end

    % Mixing without hired-in labour in Year 2
    if case_w2(i) == 2
        if case_bl(i) == 6
            psi_h(i) = w_bl(i)/h_bl(i);
            psi_l(i) = (w_w2(i) - psi_h(i)*h_w2(i))^2/(psi_h(i)*l_w2(i)^2);
            A(i) = (psi_l(i)*l_w2(i) + sqrt(psi_h(i)*psi_l(i))*h_w2(i))/(fk_w2(i)*beta*l_w2(i)^(beta-1));
        end
    end
end

% Replace missing psi_h with maximum

psi_h(psi_h==0) = max(psi_h);

% Rest of mixing without hired-in labour in Year 2
for i = 1:N
    if case_w2(i) == 2
        if case_bl(i) == 7
            psi_l(i) = (w_w2(i) - psi_h(i)*h_w2(i))^2/(psi_h(i)*l_w2(i)^2);
            A(i) = (psi_l(i)*l_w2(i) + sqrt(psi_h(i)*psi_l(i))*h_w2(i))/(fk_w2(i)*beta*l_w2(i)^(beta-1));
        end
    end
end

median_psi_l_stups=median(psi_l);
load('intermediate/Calibrated_A_psil_nonstups.mat','median_psi_l_rich');
ratio_psi_l = median_psi_l_rich/median_psi_l_stups;

%----------------------------------------------
% Estimate h*, l*, h' at rich ss level of ALL ASSETS MINUS LAND
%----------------------------------------------

k=43701;   % Upper mode of distribution of assets minus land (OB email 3/10/19)
fk=a*(k^2) + b*k;

psi_l_sim = ratio_psi_l.*psi_l;
psi_h_sim = ratio_psi_l.*psi_h;

l_opt = zeros(N,1);
h_opt = zeros(N,1);
hiredin_opt = zeros(N,1);
profit_opt = zeros(N,1);

for i=1:N
    profit = @(x) -A(i)*fk*(x(1) + x(3))^beta - w_bl(i)*x(2) + phi*w_bl(i)*x(3) + 0.5*(sqrt(psi_l_sim(i))*x(1) + sqrt(psi_h_sim(i))*x(2))^2;
    [argmax, max] = fmincon(profit, [0,0,0], [1,1,0], R, [], [], [0,0,0], [R,H(i),Hire_in_max]);
    l_opt(i) = argmax(1);
    h_opt(i) = argmax(2);
    hiredin_opt(i) = argmax(3);
    profit_opt(i) = -max;
end

clear max argmax

% Check if it would be optimal to liquidate capital:

for i=1:N
    profitalt = @(x) -rho*k -w_bl(i)*x + 0.5*psi_h_sim(i)*x^2;
    [argmax, max] = fmincon(profitalt, 0, [], [], [], [], 0, H(i));
    if -max > profit_opt(i)
        l_opt(i) = 0;
        hiredin_opt(i) = 0;
        h_opt(i) = argmax;
        profit_opt(i) = -max;
    end
end

clear max argmax

% Trim :
l_opt(l_opt < 1) = 0;
h_opt(h_opt < 1) = 0;
hiredin_opt(hiredin_opt < 1) = 0;

% Optimal cases:
case_opt = zeros(N,1);

for i=1:N
    if l_opt(i) > 0 && h_opt(i) > 0 && hiredin_opt(i) > 0
        case_opt(i) = 1;
    elseif l_opt(i) > 0 && h_opt(i) > 0 && hiredin_opt(i) == 0
        case_opt(i) = 2;
    elseif l_opt(i) > 0 && h_opt(i) == 0 && hiredin_opt(i) > 0
        case_opt(i) = 3;
    elseif l_opt(i) > 0 && h_opt(i) == 0 && hiredin_opt(i) == 0
        case_opt(i) = 4;
    elseif l_opt(i) == 0 && h_opt(i) > 0 && hiredin_opt(i) > 0
        case_opt(i) = 5;
    elseif l_opt(i) == 0 && h_opt(i) > 0 && hiredin_opt(i) == 0
        case_opt(i) = 6;
    elseif l_opt(i) == 0 && h_opt(i) == 0
        case_opt(i) = 7;
    end
end

%----------------------------------------------
% Estimate h*, l*, h' at baseline capital level
%----------------------------------------------

l_opt_bl = zeros(N,1);
h_opt_bl = zeros(N,1);
hiredin_opt_bl = zeros(N,1);
profit_opt_bl = zeros(N,1);

for i=1:N
    profit = @(x) -A(i)*fk_bl(i)*(x(1) + x(3))^beta - w_bl(i)*x(2) + phi*w_bl(i)*x(3) + 0.5*(sqrt(psi_l(i))*x(1) + sqrt(psi_h(i))*x(2))^2;
    [argmax, max] = fmincon(profit, [0,0,0], [1,1,0], R, [], [], [0,0,0], [R,H(i),Hire_in_max]);
    l_opt_bl(i) = argmax(1);
    h_opt_bl(i) = argmax(2);
    hiredin_opt_bl(i) = argmax(3);
    profit_opt_bl(i) = -max;
end

clear max argmax

% Check if it would be optimal to liquidate capital:

for i=1:N
    profitalt = @(x) -rho*k_bl_pre_trans(i) -w_bl(i)*x + 0.5*psi_h(i)*x^2;
    [argmax, max] = fmincon(profitalt, 0, [], [], [], [], 0, H(i));
    if -max > profit_opt_bl(i)
        l_opt_bl(i) = 0;
        hiredin_opt_bl(i) = 0;
        h_opt_bl(i) = argmax;
        profit_opt_bl(i) = -max;
    end
end

clear max argmax

% Trim :
l_opt_bl(l_opt_bl < 1) = 0;
h_opt_bl(h_opt_bl < 1) = 0;
hiredin_opt_bl(hiredin_opt_bl < 1) = 0;

% Optimal cases:
case_opt_bl = zeros(N,1);

for i=1:N
    if l_opt_bl(i) > 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) > 0
        case_opt_bl(i) = 1;
    elseif l_opt_bl(i) > 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) == 0
        case_opt_bl(i) = 2;
    elseif l_opt_bl(i) > 0 && h_opt_bl(i) == 0 && hiredin_opt_bl(i) > 0
        case_opt_bl(i) = 3;
    elseif l_opt_bl(i) > 0 && h_opt_bl(i) == 0 && hiredin_opt_bl(i) == 0
        case_opt_bl(i) = 4;
    elseif l_opt_bl(i) == 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) > 0
        case_opt_bl(i) = 5;
    elseif l_opt_bl(i) == 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) == 0
        case_opt_bl(i) = 6;
    elseif l_opt_bl(i) == 0 && h_opt_bl(i) == 0
        case_opt_bl(i) = 7;
    end
end

%----------------------------------------------
% Estimate value of misallocation
%----------------------------------------------

misallocation = zeros(N,1);

for i=1:N
    if case_opt(i) == 5 || case_opt(i) == 6   % Misallocation is zero if case 5 or 6 dominates at high SS
        misallocation(i) = 0;
    else
        misallocation(i) = profit_opt(i) - profit_opt_bl(i);
    end
end

misallocation(misallocation < 0) = 0;

% Plot individual-level misallocation

h1 = histogram(misallocation);

[f,xi] = ksdensity(misallocation);
figure
plot(xi,f);

misallocation_5_pcile = prctile(misallocation,5);
misallocation_95_pcile = prctile(misallocation,95);
misallocation_trimmed5=misallocation;
temp1=misallocation_trimmed5<=misallocation_5_pcile;
misallocation_trimmed5(temp1,:)=[];
temp2=misallocation_trimmed5>=misallocation_95_pcile;
misallocation_trimmed5(temp2,:)=[];

[f,xi] = ksdensity(misallocation_trimmed5);
figure
plot(xi,f);

total_misallocation_untopcoded=sum(misallocation);

% Top-code top 5% outliers
misallocation_untopcoded = misallocation;
misallocation(misallocation>misallocation_95_pcile) = misallocation_95_pcile;
total_misallocation=sum(misallocation);

% Calculate total number of individuals misallocated
num_misallocated = sum(misallocation(:) ~= 0);


%----------------------------------------------
% Calculate cost for all HHs to reach unstable ss, i.e. overcome poverty trap (at baseline pre-transfer)
%----------------------------------------------

min_transfer = zeros(N,1);

for i=1:N
    if case_opt(i) == 5 || case_opt(i) == 6   % Optimal transfer is zero if case 5 or 6 dominates at high SS
        min_transfer(i) = 0;
    else
        min_transfer(i) = uns_threshold - k_bl_pre_trans(i);
    end
end

min_transfer(min_transfer < 0) = 0;

total_transfer=sum(min_transfer);

%----------------------------------------------
% Export results to Excel
%----------------------------------------------

T=table(hhid,case_bl,case_opt,misallocation,min_transfer);
T(1:2,:);
filename='output/Structural_Results_at_ss_assets_minus_land_of_rich.xlsx';
writetable(T,filename);

save('output/Structural_Results_at_ss_assets_minus_land_of_rich')
